sumlog - #48

@devmotion I guess we also need some ChainRules methods...

@devmotion I guess we also need some ChainRules methods...

We could, but it's not necessary to do in this PR IMO.

Apart from the element type (see discussion above), I think the main problem left is that I assume the code is problematic for GPU arrays. Other array implementations in LohExpFunctions are written in a GPU-friendly way and should work with all array types.

I think the main problem left is that I assume the code is problematic for GPU arrays. Other array implementations in LohExpFunctions are written in a GPU-friendly way and should work with all array types.

Do you see a nice way of doing this?

I see two other potential things to add:

1. Support for Tuples and NamedTuples
2. Support for calling sumlog(f, x)

Both less critical, so we can come back to them

One of these updates along the way killed all of the performance. Are you seeing this too? Need to backtrack a bit I guess, and maybe split the preprocessing into a separate function

Got it

It got faster! Check it out:

julia> x = rand(1000);

julia> @btime sum(log, $x)
  6.362 μs (0 allocations: 0 bytes)
-1027.6

julia> @btime sumlog($x)
  986.857 ns (0 allocations: 0 bytes)
-1027.6

The function has to be added to the docs for the tests to pass.

I don't understand what you mean by this. Docstrings are usually added automatically, what's left to do?

I changed it to

function sumlog(x)
    T = float(eltype(x))
    _sumlog(T, values(x))
end

There's no need to restrict the type of x, and this allows it to be a Tuple or NamedTuple. For NamedTuples, calling values(x) makes it much faster, and this doesn't affect other types.

What's left to do?

You have to add the function to docs/src/index.md.

This also points to another application - maybe you really want to just compute a product, but you'd like to avoid underflow and overflow. So for example @cjdoris's LogarithmicNumbers.jl would seem to benefit from adding

Base.prod(::Type{ULogarithmic}, x) = exp(ULogarithmic, logprod(x))

avoid Base._return_type

Yes, this is part of the code complexity goal but will also improve stability of LogExpFunctions. All such internal functions and "hacks" should be removed from the PR, in particular since it seems they can be avoided easily. Even standard libraries such as Statistics don't use _return_type to handle empty iterators, see eg https://github.com/JuliaLang/Statistics.jl/blob/cdd95fea3ce7bf31c68e01412548688fbd505903/src/Statistics.jl#L204 and https://github.com/JuliaLang/Statistics.jl/blob/cdd95fea3ce7bf31c68e01412548688fbd505903/src/Statistics.jl#L170.

@cscherrer: Regarding stepping back and agreement: I always think in terms of costs and benefits (code complexity and maintainability vs how useful the code is), and personally I would just go to logs as soon as possible, even at a slight performance cost. But if you really need this, I am fine with including it.

Regarding the goals:

1. I would rename to logprod, conveys the underlying algorithm and precision trade-offs better,

2. I think that dims and foo(f, x) are unnecessary, and it is silly that each function replicates the boilerplate for this, given that Julia has much nicer mechanisms now for these, but I understand that whenever we leave these out someone will complain

3. I would be happy with a robust and reasonably accurate logprod that is approximately sum(log, ...), with the understanding that someones one is better than the other. All algorithms have trade-offs and that's fine. Maybe we should document them though.

I've pushed another version. This time

1. It's logprod instead of sumlog
2. sumlog is still there for now, for easy comparison
3. I remembered about frexp. @mcabbott this behaves better for subnormals
4. I dropped dims, etc. I think we have a better understanding now of each other's priorities. I'm still in favor of more functionality, but we can start simple and get tests etc going. That will make it easier to weigh any drawbacks of adding more functionality.
5. I think this is a good candidate to go in Base. If we make it logabsprod it can be a big help to speed up logabsdet:

julia> x = LowerTriangular(randn(1000,1000));

julia> using LinearAlgebra

julia> using BenchmarkTools

julia> @btime logabsdet($x)
  8.687 μs (0 allocations: 0 bytes)
(-631.836, -1.0)

julia> d = diag(x);

julia> @btime logabsprod($d)
  1.202 μs (0 allocations: 0 bytes)
(-631.836, 1.0)

I cheated here a little, since we don't have (that I know of) a lazy diag in Base.

logprod is a neat idea to avoid checks.

frexp is also clearly what we were looking for. Does this have any effect on speed?

I am lost in all the noise on minor details, but checking sig > floatmax(typeof(sig)) / 2 is now the wrong thing, as it will overflow towards zero.

logprod should really have a case which tries log(prod()) first on small enough arrays, as this is much faster. (And it should advertise itself as being less prone to overflow than log(prod, rather than as being faster.)

logprod is a neat idea to avoid checks.

Thanks! I think the name is more natural too :)

frexp is also clearly what we were looking for. Does this have any effect on speed?

I had assumed if anything it might be more efficient, but maybe not. Here's a quick check:

julia> @benchmark frexp(x) setup=(x=100 * rand())
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     2.063 ns (0.00% GC)
  median time:      2.084 ns (0.00% GC)
  mean time:        2.085 ns (0.00% GC)
  maximum time:     6.112 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000

julia> @benchmark (significand(x), exponent(x)) setup=(x=100 * rand())
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     1.613 ns (0.00% GC)
  median time:      1.633 ns (0.00% GC)
  mean time:        1.640 ns (0.00% GC)
  maximum time:     10.330 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000

May not be a real effect since this is sub-nanosecond.

I am lost in all the noise on minor details, but checking sig > floatmax(typeof(sig)) / 2 is now the wrong thing, as it will overflow towards zero.

I had missed this, but @devmotion caught it too. Very surprising it would be so different.

I'm also wondering... Was this wrong in the first place? A factor of 2 should work if we assume it's sequential. But some implementations might exploit the associativity. In this case I think we need sqrt, or better, to look at exponent(sig). But with that we end up in the renormalization branch more often, so we're only about twice as fast as sum(log, x)

logprod should really have a case which tries log(prod()) first on small enough arrays, as this is much faster. (And it should advertise itself as being less prone to overflow than log(prod, rather than as being faster.)

Both great points.

do we want to merge this?

I lost track of the various changes, but I am fine with merging. We can always micro-optimize things later; this is already more efficient than sum(log, x).

Same. I expected this would be relatively simple, and was surprised by the number of cases to worry about. I agree merging and then handling various cases as they come up is reasonable. It looks like the tests need some update though, they're currently failing because they use sumlog (now undefined) instead of the updated name logprod.

I don't remember the details either. It seems there are unaddressed comments and the tests don't pass yet, so I guess the PR needs some updates and another round of review before merging.

Friendly ping: this is a great contribution and it would be unfortunate to leave it dormant. @cscherrer, if you have the time to address the May 10 review comments by @devmotion and fix CI, I would be happy to merge.

Closing because of no activity.

That said, this is useful functionality, so if there is in interest in reviving the feature, just please open a new PR.